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We present a theory of dark soliton dynamics in trapped quasi-one-dimensional Bose-Einstein 
condensates, which is based on the local density approximation. The approach is applicable for 
arbitrary polynomial nonlinearities of the mean-field equation governing the system as well as to 
arbitrary polynomial traps. In particular, we derive a general formula for the frequency of the 
soliton oscillations in confining potentials. A special attention is dedicated to the study of the 
soliton dynamics in adiabatically varying traps. It is shown that the dependence of the amplitude 
of oscillations vs the trap frequency (strength) is given by the scaling law Xo oc uj ' where the 
exponent 7 depends on the type of the two-body interactions, on the exponent of the polynomial 
confining potential, on the density of the condensate and on the initial soliton velocity. Analytical 
results obtained within the framework of the local density approximation are compared with the 
direct numerical simulations of the dynamics, showing remarkable match. Various limiting cases are 
addressed. In particular for the slow solitons we computed a general formula for the efi'ective mass 
and for the frequency of oscillations. 
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I. INTRODUCTION 

One of the main properties of solitons, making them to 
be of special interest for physical applications, is preserv- 
ing their localized shapes during evolution and mutual 
interactions 0. Due to this robustness solitons can be 
regarded as quasiparticles and systems possessing large 
number of such excitations can be described in terms of 
the distribution function governed by the kinetic equa- 
tion 0. 

In the mean-field theory 1^ description of the quasi- 
one-dimension homogeneous Bose gas is reduced to the 
exactly integrable nonlinear Shrodinger (NLS) [or one- 
dimensional (ID) Gross-Pitaevskii (GP)] equation, and 
therefore solitons are expected to play a prominent 
role in the dynamical and statistical properties of low- 
dimensional condensates. When interatomic interactions 
are repulsive, the GP equation possesses dark (or grey) 
soliton solutions 0, 0, fs] . Existence of the dark solitons 
was confirmed by a number of recent experiments with 
BEC's confined by elongated traps 0. 

In practice, condensates appear to be never homoge- 
neous, and therefore effect of external potentials on the 
dark-soliton dynamics is a subject of special interest (see 
e-g- IE Hi and references therein). An inhomo- 
geneity of a system by itself does not invalid possibili- 
ties of description of solitons as quasiparticles (in some 
approximation, of course). In particular, one can ex- 
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plore Hamiltonian approach to an effective particle with 
one degree of freedom, instead of dealing with the orig- 
inal equation for the macroscopic wave function, which 
is a system with infinite degrees of freedom. Moreover, 
one can extend the respective description on the gas of 
solitons, which now will be described by a distribution 
function governed either by Fokker-Planck equation (for 
the case where a soliton bearing system interacts with a 
thermal bath, see e.g. 'iS*!) or by a kinetic equation with 
respective collision iirtcgral, as this is shown in for 
the case of interaction of solitons with a noncondensed 
atoms. 

A quasiparticle description of dark solitons can be ob- 
tained from the perturbation theory in adiabatic approx- 
imation Pi (sometimes called the collective variable ap- 
proach). At the same time, as was shown in [Tol |. a con- 
cept of a quasiparticle naturally emerges from the Lan- 
dau theory of superfluidity and can be justified on the 
basis of the mean field theory within the framework of 
the local density approximation. It turns out, that a dark 
soliton moves in an external potential without deforma- 
tion of its density profile as a particle of mass 2m. The 
local density approximation is rather general, allowing 
direct extension to other nonlinear equations, related to 
the BEG dynamics, as well as to various (non-parabolic) 
types of the trap potential. Building up such a general- 
ized theory is the main goal of the present paper. 

In real experimental conditions the external trap po- 
tential can depend not only on coordinate, but also on 
time. That is why the second aim of the present work is 
the description of the effect of adiabatic time-dependence 
of the external parameters on the dark soliton motion. 

The paper is organized as follows. We start with the 
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dynamics of a dark soliton in an adiabatically chang- 
ing parabohc trap (Sec. In Sec. 11111 we develop our 
Hamiltonian theory for solitons described by generalized 
polynomial NLS equations and show how such approach 
is related to the mean field approximation. In Sec. IIVI 
we consider in detail examples of dark soliton dynamics, 
which include the cases of non-parabolic trap and models 
with higher nonlinearity. The consideration is provided 
within the framework of the local density approximation 
and is verified by direct numerical simulations of the dark 
soliton dynamics. In this section we also show how one 
can modify the perturbation theory for dark solitons to 
take into account adiabatic change of the trap frequency 
('Sec. lVI3|l and make comments on the dynamics of small 
amplitude dark solitons (Sec. IV C|l . The outcomes are 
summarized in the Conclusion and technical details of 
some calculations are given in the Appendices. 



II. DARK SOLITON IN A TIME-DEPENDED 
PARABOLIC TRAP 

Let us start with the dynamics of a dark soliton de- 
scribed by the GP equation 



2m 



^ + g\^\^^ - fi^. (1) 



Here g = 2h^as/{ma\), as is s-wave scattering length 
and a_L the transverse linear oscillator length, which de- 
scribes the BEC in an elongated trap at low densities 3j 
(see also for the details of derivation by means of the 
multiple scale expansion method). 

It has been shown in Ref. JJ]] (see also the details below 
Sec. IIII|I that the dark soliton dynamics in a parabolic 
trap can be successfully described within the framework 
of the local density approximation. This means that, 
in spite of the presence of the trap, one starts with the 
solution of the ID homogeneous (i.e. when lu = 0) GP 
equation |^ (see also 0], §5.5): 



tanh 



x~X{t) 



,(2) 



where X{t) = vt, v is the velocity of the soliton, uq is the 
unperturbed linear density, c 
sound and £ = h/ {m^/TF' 



griQ/m is the speed of 



v"^) is the width of the soliton. 
Then the influence of the trap is accounted by considering 
a general function X{t) which dependence on time is to 
be obtained. 

The energy of the system can be defined as 



E = 



2m 2 ^ ' 



hen, (1 - ^ 



dx 

2\ 3/2 



(3) 



and for the dark soliton solution (0) can be rewritten in 
a form of the conservation law 



(X) 



= {QE) 



2/3 



(4) 



where Q = Zg/{Ahm). The introduced dependence c = 
c{X) is the key point of the local density approximation: 
the sound velocity is substituted by its local value com- 
puted in the point where the center of the soliton is lo- 
cated. In the Thomas-Fermi (TF) approximation, when 

the atomic density is given by n{x) = ^ ^/i — "^"^ ^ ^ , 

one has 



c'{X) = -n{X) 
m 



loj^X^ 



with Co = \J /i/m. Substituting v — dX/dt in E 
energy conservation can be rewritten as follows [IC 
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dX 
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Here we introduced the effective mass of the soliton con- 
sidered as a quasiparticle 



2m, 



the frequency of the soliton oscillations = lo/ \/2 
llOl IllL [1^ and the effective soliton energy 



(7) 



E^ = 



i5E) 



2/3 



(8) 



which, altogether with E, is a constant of motion. The 
amplitude of oscillations governed by (jBJ is 



Xq — \j2E^j 



(9) 



One of the characteristic features of the introduced 
quasiparticles is that their dynamics is determined not 
only by their local properties (velocity and amplitude) 
but also by the environment, i.e. by the unperturbed 
density. As a result any change of the trap characteris- 
tics (say, trap frequency or geometry) will affect solitons 
not only by changing the domain of their motion but also 
through the change of the density. It turns out that the 
local density approximation is a suitable framework for 
description of mentioned phenomena in the case when 
time variation of the parameters of the system is slow 
enough. 

According to a general law of the Hamiltonian mechan- 
ics, the adiabatic invariant 
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stays constant 14]. Time dependence of the amplitude 
of oscillation can be defined from this condition. The 
canonical momentum, which enters in ()10|l . can be com- 
puted explicitly using the formula 



P = 



" dE_(h}_ 

dv V 



(11) 
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what gives 




arcsm 



(12) 



It turns out, however, for calculation of the adiabatic 
invariant it is more convenient to use the general equation 
between / and the frequency of oscillations: 



dl I f dX , 



(13) 



Taking into account that tUg = w/ V2 does not depend 
on E and using an obvious boundary condition / = at 
w — > 0, we easily find a simple equation 



(14) 



It is not difficult to show (see, for example, Eq. (17.10) 
in Q), that in the TF approximation one has 

Cq — gn{0)/m(x oj"^^"^, i.e. Cq cx w , (15) 
so the second term on r.h.s. of H14|) is constant. 



Thus preserving the adiabatic integral in an adiabatic 
process implies preserving E/uj, what in the case of 
slowly varying frequency implies E oc u. Taking again 
into account that according (|15|l in the TF approxima- 
tion Co oc w^/'^, one deduces from ^ that cx u;'^/^. 
Finally, the scaling law for the amplitude of oscillations, 
defined by reads 



Xo cx 



(16) 



It is worth to underline that this law is different 
that one for a conventional harmonic oscillator, where 
Xq (X tj"^/^, even though the motion of the soliton is 
pure harmonic. The point is that in our case the ratio 
E/lu, but not E^/uj, is preserved. 

An important feature of the soliton dynamics is that in 
the case at hand the soliton frequency does not depend on 
the energy. Hence the frequency of the soliton oscillations 
does not depend on the amplitude of the soliton, what 
corroborates with the analysis of the oscillations of the 
small- amplitude solitons [see (|^ below and subsequent 
discussion] as well with the earlier studies 0, ■ 

We have checked the obtained predictions, made on 
the basis of the local density approximation, numerically. 
The typical results are presented in Fig. ^ 



The local density approximation essentially uses that 
the background of the condensate is static, i.e. that the 
dark soliton motion does not excite the motion of the 
whole condensate. In practice, due to finiteness of the 
system, such a supposition strictly speaking does not 
hold, and the whole condensate also undergoes oscilla- 
tions with the frequency of the condensate iv, what fol- 
lows directly from the Ehrenfest theorem. The difference 
of the frequencies of the condensate and of the dark soli- 
ton, i.e. between to and tOs, results in the beating of the 
dark soliton 0, which are clearly observable in Fig. ^ 
Respectively, one can identify the two slopes correspond- 
ing to the maxima and to the minima of the soliton am- 
plitudes. We will use the subindexes and "— " for the 
respective quantities. In other words, each of the results 
presented in panels (a) - (c) and (e) - (g) is character- 
ized by the two scaling laws: Xq,± = X-i-o;"'''* shown 
explicitly in Fig. ^ (a),(e). The exponents 7± are differ- 
ent (although their difference is relatively small), what 
requires a definition of some averaged exponent 7 which 
could be compared with the theoretical predictions. We 
obtain such exponent numerically from the dynamics of 
the averaged amplitude, i.e. using the formula 

Xq 1 f X+ X_\ 

— = 77 — + • 17 



the figures the law of the change of the amplitude of soli- 
ton oscillations stays close to the predicted law 7 = 2/3 
for relatively slow solitons and relatively large densities. 
Meantime deviations are clearly seen in Fig.n(c) and (e). 
In the last case the exponent 7 is essentially less than the 
predicted in our analytical consideration. It turns out, 
however that the mentioned deviation from 2/3 law is ob- 
served for small densities. This is natural from point of 
view of the theory. Indeed, our consideration was based 
on the TF approximation for the atomic density, when 
no oc uj^^^. This approximation fails at low densities, and 
must be substituted by the Gaussian distribution, where 
no cx uj^/^. Then by repeating the above our arguments 
for the Gaussian distribution, instead if the TF one, one 
finds 

Xqcxw-i/^ (18) 

i.e. the law of the dependence of the amplitude of oscilla- 
tion of the conventional linear oscillator on the frequency, 
what corroborates with the numerical findings. 



III. GENERAL APPROACH 
A. Generalized equation. 



Summary of the results for the averaged exponent 7 are 
presented in panels (d) and (h). As one can see from 



The theory developed in the previous section can be 
generalized for NLS equation with arbitrary power-law 



4 




FIG. 1: Time dependence of the soliton coordinate on the adiabatically changing frequency, modeled by the function uj{t) = 
(1 + O.OOlt) ■ 0.1, in logarithmic coordinates: panels (a)-(c) and (e)-(g). Straight solid lines visualize the law Xq = Xo^ '• 
Dashed lines in panels (a) and (e) show the exponents 7± [see H17|l ] while in the rest of panels only the average 7 is shown. 
In panels (a)-(c) the parameters are no = 1 and v = 0.1; 0.5 and 0.8, correspondingly. The respective matching parameter 
is Xo = 0.385; 1.93; 2.6. In panels (e)-(g) the parameters are v = 0.1, no = 0.3; 0.6, and 1, correspondingly. The respective 
matching parameter is Xq — 0.622; 0.445; 0.385. In panel (d) we plot dependence of 7 on the initial soliton velocity v for the 
case no = 1. In panel (h) we depict the dependence of 7 on the initial density no corresponding to the case of the initial velocity 
V — 0.1. In the numerical calculations we take h — 1, m — 1 and g = 1. 



nonlinearity and non-parabolic potential. More specifi- 
cally, in the present section we consider the equation 

z;i*t = -^^',, + ;7(x)^' + g|*|2"*-^* (19) 
Zm 

where a is a positive integer and <? > 0, which describes 
interacting particles of mass m in an external poten- 
tial U{x). The exponent a characterizes the effective 
inter-particle interactions. In particular when a = 1 and 
g = 2h'^as/ {ma^) one recovers the GP equation Q con- 
sidered in the previous section. 

The chemical potential ^ introduced in (|19|l is deter- 
mined by the link valid for a homogeneous condensate: 
fi = gn°'. Thus the sound speed c connected to the chem- 
ical potential by the relation mc^ — ndfi/dn can be ex- 
pressed as follows 



ag 



(20) 



There are several reasons to consider more general 
equation H19|l . First of all, equation being com- 
pletely integrable, possesses very specific soliton proper- 
ties. It is interesting to investigate the soliton dynamics 



in a more general situation. The case of a = 2 is par- 
ticularly important, because corresponding equation can 
be used in different physical problem. Such a situation 
can take place near the Feshbach resonance. In this case 
the s-wave scattering length depends on magnetic field 
as flg — ag + A/{B — Bq) where ag is the background 
value of the scattering length, and Bq and A are the lo- 
cation and width of the resonance. If magnetic field is 
equal to Be = Bq — A/og, the scattering length turns to 
zero and the dominant interaction among atoms is due 
to three-body effects. 

Indeed, in the higher approximations of the Bogoliubov 
theory expansion of the chemical potential of an uniform 
gas with respect to density n has form 



fj, — agTi 



bi + b2inaiy/^ + bsinal) In 

nai 



52n',(21) 



where hi = AttIt^ /m and other coefficients b can be calcu- 
lated (see 0, §4.2). Coefficient 32 depends on three- 
body interactions and cannot be calculated explicitly. 
However, it stays finite for B = B^,, while three first 



terms disappear, giving /i = uM- Correspondingly 
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the non-linear term in the mean-field equation has form 
52 I ^' 1^ The sign of g2 cannot be defined from general 
considerations. We assume that 32 > 0. After averaging 
with respect to the transverse motion we obtain H19|) with 
a = 2 and g = g2 / [ii^^ a\_) . 

Another physical system where the equation of the 
state with a = 2 is valid, is an ID Bose gas in so- 
called Tonks-Girardeau (TG) limit of inpenetrable parti- 
cles. This limit can be achieved for a gas of small density. 
It has been shown by Girardeau that there exists an ex- 
act mapping between states of this system and an ideal 
ID Fermi gas. In particular in this case one has fi = gn? 
with g = ?i^7r^/(2m). It has also been rigorously shown 
that one can find density distribution of a such gas in a 
ID trap by minimization of the energy functional 




On the basis of these considerations authors of Ref. 
[18] suggested to use equation (19) for dynamics of the 
TG gas. However, the hydrodynamic-like equation (19) 
can not give a satisfactory description of dynamics of 
an ideal Fermi gas. Nevertheless it can be useful for 
a Bose gas near the TG limit, where equation of state 
approximately follows the a = 2 law, but dynamic is still 
not an ideal gas- type. 

The case a = 2 is often referred also as a quintic nonlin- 
ear Schrodinger (QNLS) equation. For the sake of brevity 
in what follows use this terminology. 

We mention that other polynomial models are also con- 
sidered in literature [l^. 

B. Soliton in the generalized equation. 

Let us consider now a condensate in the absence of 
external field, U{x) = 0. Eq. (|19|l takes the form 

i?i*t = -^«',,+5|*P"*-^* (23) 
Zm 

and is subject to the finite density boundary conditions: 
lim *(2;,t) = V^e=^'^ (24) 

X — ^±00 

where the constant 9 can be considered without restric- 
tion of generality in the interval [0,7r/2]: 9 G [0, 7r/2]. 
Then dark solitons, ^'s(a;, t), will be associated with trav- 
eling wave solutions, characterized by the following de- 
pendence of the density on the spatial coordinate and 
time: 

\^,ix,t)\^=r^^ix-vt), (25) 

where v is the soliton velocity. Below such solutions will 
also be referred to as unperturbed. 

The energy of the soliton solution can be defined as 

I £{x)dx (26) 



where the energy density £{x) is given by 

^5«o"(l*l'-"o) . (27) 

The energy is an integral of motion. Hence taking into 
account that the dark soliton depends on the two param- 
eters (no,i;) and connecting the mean density with the 
speed of sound by (pn|l . one concludes that the energy of 
the dark soliton is a function of c and v: 

E^E{c,v). (28) 

C. The local density approximation. 

Consider now propagation of a dark soliton in a con- 
densate with the density, varying due to the external trap 
potential: n = n{x) with n(0) = no (for the sake of defi- 
niteness the trap potential will be assumed having min- 
imum at a; = 0: J7(0) = 0). In particular, in the TF 
approximation the function n{x) is given by 

n{x) = UTFix) = - Uix)]"'' . (29) 

This formula determines the dependence of the sound 
velocity on the spatial coordinate [c.f. ©]: 

c^ix) = cl-—U{x), (30) 
m 

where cq is expressed through tiq by the link l|2()|l . 

Now we define the local density approximation as an 
assumption that the conservation law (|28|) is valid for 
a soliton in the inhomogeneous condensate, i. e., that 
c can be changed to its local value c{X), where X is 
the position of the center of the soliton, computed using 
the unperturbed soliton wave function 4's(x,i). Respec- 
tively, X and V are considered as functions of time related 
by the equation dX/dt = v(t). 

Thus, in the local density approximation the equation 
of motion of the soliton is determined from H28f) : 

E{c{X),v) E ^0. (31) 

Here E is the constant energy of the soliton. 

Eq. 131|l can be viewed as an equation of motion of a 
quasiparticle, which can be associated to the dark soliton. 
Then E(c{X),v) must be associated with the Hamilto- 
nian of the quasiparticle after expressing the velocity v 
through the canonical momentum p according to the for- 
mula Hll|) . After inverting this formula, one obtains the 
Hamiltonian of the quasiparticle: 

H{p,X) = EiciX),vip,X)). (32) 

Finally, the adiabatic invariant and the frequency are 
computed according to formulas (|10|) and (|13|l , which are 
obviously valid in the general case. 
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D. Justification of tlie local density approximation. 

In the present subsection we show that equation for the 
energy of a sohton, obtained for an uniform condensate, 
is actuahy vahd also for a trapped condensate in the focal 
density approximation. Thus the trapping potential does 
not enter explicitly in the expression for the energy of a 
soliton. 

To this end we define a real-valued wavefunction of 
the background F{x) such that no(a;) — noF'^{x) is the 
density of the condensate in the absence of the soliton 
and F{x) solves the equation 

- ^F,, + gn^F^''+' + [U{x) ~ ^i]F ^ (33) 

where no = no(0) subject to the normalization conditions 
F{0) = 1 and F^{0) = 0. 

The density of the "grand canonical energy" of the in- 
homogeneous condensate can be written down as follows 

= I^I^.P + jj-n^+Hx) + [U{x) - fi]n{x) . 
Zm I + a 

(34) 



Here n{x) = Let the soliton center he a,t x = X , £ 

be a soliton width and Lq be the spatial extension of the 
condensate. Then we introduce S such that Lq S ^ £ 
and separate the integration on two domains 

E' ^ J £'dx+ j £'dx. (35) 

\x-X\>S \x-X\<5 

Next, we add to the first term an integral /|a;_x|<5 '^o'^-^ 
and, correspondingly, deduct it from the second term in 
E'. 

For the case of a dark soliton solution, which is expo- 
nentially localized around x = X, the first integral can be 
approximated (with the exponential accuracy) as follows 

j £'dx+ j £'^dxf=i j £^dx^Eo (36) 

\x-X\>S \x~X\<S °° 

where Eq is the energy of the unperturbed condensate. 

In order to compute the other two integrals we repre- 
sent 



£'{x)-£^{x) 



2m' 



1 



[n"+i(a;) - n^+\x)] - gn^{x) [n{x) - no(a;)] + 



f^F^ 
2m F 



■n{x) 



4m 



iP')xx -(37) 



As it is shown in Appendix^ the last two terms can be made as small as necessary by choosing the potential large 
enough, while in the rest of the terms related to the background x can be securely substituted by X (due to their 
smoothness in the region of the soliton motion). This leads us to the final expression for the energy of the soliton: 



E, = 



{£' - S^) dx 



l*-l' + ^[l*l'"^'-"o+'W]-5«oW \^f-noiX) }dx.{38) 



2m 



a + l 



\x-X[<S 



\x-X\<5 



The obtained integral does not depend (in the lead- 
ing order) on the particular choice of the parameter 5. 
Then comparing the expression H38|) with (|25|) . H27|) one 
can verify that they lead to the same expression for the 
soliton energy, where the only substitution no by no (A") 
must be made. 



IV. EXAMPLES OF LANDAU DYNAMICS OF 
DARK SOLITONS 

In the present section we consider two examples rele- 
vant in different ways to the BEC dynamics in low di- 
mensions. 



A. Dark soliton of the GP equation in a 
polynomial trap. 

1. General approach. 

Let US now turn to the case where the "polynomial" 
trap 

with r being a positive integer, r = 1, 2, and u being 
a function slowly depending on time: uj — uj(t). If r = 1, 
then U{x) is transformed in the conventional parabolic 
trap considered in Sec. ^ Then uj is the trap frequency. 
For this reason and for the sake of brevity of notations in 
what follows uj is referred to as a frequency independently 
on the value of r. 

The question we are interested in is the dependence 



7 



of the amplitude of the sohton oscillations on the fre- 
quency, subject to the adiabatic change of the last one. 
The explicit form of p, given by p2|l , allows one to solve 
the problem analytically in a general case, i.e. for the 
arbitrary integer r. 

Now the link between the velocity and the coordinate 
Q reads 



^2 + 1^2r^2r ^ ^2 



(40) 



[c* was defined in ||HJ)] and expression for the amplitude 
of the oscillations of the soliton, Xq is given by: 



Xn 



2l/(2r)^l/'- 



(41) 



where 



Rr — 



21/(2,-) 



dx 



1 vT 



„2r 



(46) 



The obtained relations, as well as constancy of the total 
number of particles N and of the adiabatic invariant / 
subject to slow change of the frequency readily allow one 
to get the scaling relations [they follow from H42f) and 
H44|l . respectively]: 



no 



(47) 



Next one can compute the following quantities: 
- The normalization condition 



= / n{x)dx 

-Xtf 



2r(2no)i+i/(2'^) / g \ 1/(2^) 1 



2r + 1 



(-) 

Vm/ 



(42) 



where N is the total number of atoms and we introduced 
the TF radius 



XTF 

The adiabatic invariant 



(43) 



(44) 



where the constant Gr is defined in (|C2(I and the details 
of calculations are presented in Appendix lUl 
- The frequency of the soliton [using 



(45) 



Finally, taking into account the link (|41|l we arrive at the 
general scaling relation determining the dependence of 
the amplitude of the soliton oscillations on the frequency 



Xq OC LU 



2r 
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1 + 2r 

2. GP dark soliton m an x'^-trap. 



(48) 



Let us consider in more details dynamics of a soliton 
in a non-parabolic trap with the potential energy 



U{x) = ^uj^x^ 



(49) 



(i.e. the case r = 2). Now R2 = tt2^/'^K{1/V2) « 0.847, 
K(-) being the complete elliptic integral of the first kind, 
and the frequency of soliton oscillations depends on the 
energy of the condensate [see (|40|) and 145() ]. The expo- 
nent defined by (gHI> is 7 = 0.8. 



The numerical study of the soliton dynamics in a quar- 
tic trap are presented in Fig. |2J While the predicted 
exponential law uj~^'^ is now also obtained with reason- 
able accuracy, there are several features which distinguish 
the present case from the case shown in Fig. ^ First, 
one does not observe beatings (while they are well pro- 
nounced in the case of a parabolic trap) . This fact can be 
explained by the absence of the unique frequency of the 
background oscillations: in the case at hand the Ehren- 
fest theorem does not result in a coupled equation for the 
averaged coordinate of the center of mass of the conden- 
sate. Second the dependencies of the exponent 7 of the 
soliton velocity and on the density appear to be decreas- 
ing functions, as it is shown in Figs. [21(d) and (h). 



B. Dark soliton in the QNLS limit. 

As the next example we consider the equation 



ih^t = -^"^xx + -muj^x^^ + - • (50) 

Now a = 2 and r = 1. Although general approach, simi- 
lar to one developed in the preceding section is also avail- 
able in the case at hand, it becomes rather involved and 
cumbersome. That is why, here we consider the phys- 
ically relevant case of the parabolic potential which re- 
veals the main physical features of highly nonlinear mod- 
els. 
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0.2 0.4 0.6 0.8 ■ 0.4 0.6 0.8 1 

FIG. 2: Dependence of the soliton coordinate on the frequency in the logarithmic scale. Adiabatic change of the frequency is 
modeled by the law uj{t) = (1 + O.OOlt) • 0.1. The dashed lines visualize the law Xq — Xquj In panes (a), (b) and (c) the 
parameters are no = 1 and v — 0.1; 0.3, and 0.6, respectively. The matching parameter is Xo = 0.49; 0.98; 1.52. In panels (e), 
(f) and (g) the parameters are v = 0.1, and n — 0.4; 0.7, and 1, respectively. The matching parameter is Xo = 0.38; 0.445; 0.49. 
In panel (d) we show the dependences of the exponent 7 vs soliton velocity v at the density no = 1. In panel (h) we show the 
dependences of the exponent 7 vs density no at u = 0.1. In the numerical calculations we used h = 1, m = 1 and g — 1. 



The dark soliton solution has the following form Taking into account that due to 1)29(1 now 



= v/;^:(M)e''^("'*^ (51) 
^ ^ ^ 12V5no(c2-i;2)e(-^(0)/^ 
ns(xA) — voTio 9 

c2(4 + P>-^(*))/^)'-12(c2-w2) 

(52) 

, , , /c2e(--^W)/^-2c2 + 6«2\ 
6'5(a;,i) = -arctan — (53) 



where Xify — + xq, xq is a constant, and I 
The TF distribution now acquires the form 



nTF(x) — \ u muj'^x'^ 



(54) 



and the normalization conditions defines the chemical po- 
tential = \/2rngujN /tt. 

The energy is computed from (|26|l . (|27|l to be 



i^ = ft./^^(c^-.^)ln?^±^. (55) 



g 4^/2 



2c- 



2 2 

- UJ X 



and introducing the notation 



£0 = hllQCQ 



(56) 



(57) 



we obtain 



E = ^0^(1- 



X In 



2-\/Cq — LO'^X'^ + V^\/ Cq — uj'^x'^ — v"^ 

2-\/cq — uj'^x'^ — \/3\/ Cq — uj'^x'^ — v"^ 



(58) 



u 



Respectively, the energy of the zero velocity dark soliton 
is Eo = E{v = 0) w 1.145o- 

In Fig.|3Ia) we present a typical trajectory of the QNLS 
dark soliton in a constant trap. One of the main features 
observed is that the dynamics is not strictly periodic, but 
undergoes slow modulations (see Fig. ^ . The averaged 
frequency of the dynamics shown is approximately 0.07 
(this corresponds to the relation ujg w 0.7w) while the 
frequency of the large oscillations of the period is ap- 
proximately 5 times less. It is worth pointing out that 
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the theoretical prediction for the frequency of the large 
amplitude (slow) dark solitons in the QNLS model gives 
uig ~ 0.6572ci; (see the TableQlbelow) while small ampli- 
tude solitons should oscillate with the frequency close to 
Ws « w (see the discussion in Sec. IV Cp . 



FIG. 3: (a) Dependence of the position of the center of the 
QNLS dark soliton on time for uj = 0.1, no = 1, and v = 0.1 
(b) Time dependence of the half-period T/2 substracted from 
figure (a). As before we take h — 1, m — 1 and g = 1. 
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For the next step we studied the adiabatic dynamics 
of the QNLS dark soliton in a slowly varying trap. The 
respective results are shown in Fig. 01 




FIG. 4: Dependence of the soliton coordinate on the frequency adiabatically varying according to the law Lj{t) = (1-1-0. 001t)-0. 11. 
Straight lines show the law Xq oc u> In the panels (a) and (d) we show also the laws Xo.± oc tj~^± [see llTIl l by dashed 
lines. In figures (a), (b) and (c) parameters are no = 1, and v = 0.14; 0.42, and 0.85, correspondingly. In figures (d), (e) and 
(f) parameters are v — 0.14, and no = 0.4; 0.6, and 1, correspondingly. In the numerical calculations we take h — 1, m = 1 and 
<7 = 1- 



Like in the case of the NLS dark soliton in a parabolic 
trap one can observe beating of the solution. From the 
left column [panels (a) to (c)] one detects increase of the 
frequency with increase of the initial value of the velocity, 
what is expectable in view of the above discussion. The 
right column [panels (d) to (f)] show that the frequency 
of soliton oscillations decay when the background density 



increases. This last fact is also explained in view of the 
above discussion, by the fact that increase of the local 
density subject to constant velocity v results in increase 
of the relativistic factor — w^, and thus in bigger dif- 
ference between the speed of sound and soliton velocity. 
In all the cases however one observes well pronounced 
scaling law with the exponent 7 = 0.5 for relatively large 
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range of the parameters. 

V. LIMITING CASES 
A. Small velocity solitons. 

1. General relations. 

As we have seen above, increase of the power of the 
nonhnearity (i.e. of the exponent a) makes the problem 
of the computing the frequency and dependence of the 
frequency on the ampUtude of osciUations rather comph- 
cate, not allowing one to obtain a general formula linking 
Xq and uj for arbitrary a. It turns out however, that the 
problem can be solved in the limit of small velocities: 
V <^ c. To this end we take into account that the static 
dark soliton for any a > has zero amplitude in its cen- 
ter, and hence the limit of small velocities corresponds 
to the limit of small X. Then, expanding Eq. (|31l) with 
respect to v"^ and X^, one obtains in the leading orders 



E = Eo + ^v' + ^X' 



(59) 



where the subindex "0" stands to indicate that the re- 
spective quantities are computed in the point X ~ Q and 
V = Q. This formula must be viewed as a standard ex- 
pression 



(60) 



solution, i.e. one having form H25|l and thus depending 
only on the running variable x—vt (what in dimensionless 
variables means dependence on C — 2Vt), and represent- 
ing ip = riexp{i6) one obtains (see Appendix IB|) the link 



-V 



l-rf 



(65) 



and the equation for r] (notice that according to H24|) 
the boundary conditions now are — > 1 and 0^ — s- as 
C->±oo) 



^cc + (1 - ^'") v~v 



jy3 



= 0. 



(66) 



The last equation can be integrated with respect to C 
what gives 



— (7?2"+2 - 1) + 1 - _ y2 (L^?!L (67) 

a -|- 1 rj^ 



where we designated rjQ = P = P(?/). 

Now the energy of the soliton can be rewritten in the 
form (see Appendix 



2V2 



..... , 



P2(^) 



7?- 



drj 



(68) 



where £q was introduced in (|57|l and rjm determines the 
soliton amplitude in its center and solves the equation 



for the enrgy of a harmonic oscillator having mass rUg 
and frequency ujs- Comparison of H59|l with l|60|) gives 
the expressions for the effective mass 



rus = 2— -y 
and for the frequency of oscillations 



dEp / dEp 
dX^ / dv^ 



1/2 



(61) 



(62) 



of a small amplitude dark soliton. 

Thus to compute frequency dependence of the ampli- 
tude of the soliton oscillations from Eq. (|23(l we have to 
expand the energy E{c{X),v) for small X and v. It is 
convenient to do this in dimensionless variables which we 
define as follows: 



= np ^/'vl/, 



c 



mcp 



2 

-X, 

a 



-t, (63) 



h \ a ha 
allowing one to rewrite H23|) in the dimensionless form 

i^r = -Ac + (1^1'" - (64) 

[here we used the relations H20|l ]. Also we will use the 
notation V = \f%—. Then looking for the dark soliton 



P(j]ra) = 0. 



(69) 



For a particular case of the zero-velocity dark soliton 
one has 



where 



Gq — 



2^/2 



\J a{a + 1) 



\/ 7y2a+2 „ (q, _|_ 2) 7^2 _|_ ^ 



(70) 



(Notice that in this case r\ra — 0, but P(t] = 0) = 
yja/{a + 1) 7^ 0.) Particular values of the energy for 
some relevant models are presented in Table 



2. Effective mass of a dark soliton. 

Let us consider now a soliton moving with a small ve- 
locity, V <^ 1. To execute the expansion of the energy 
we first notice that from H69() and (|67|l it follows that in 
the leading order 



-V. 



a 



(71) 
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Next, we introduce a constant 770 which satisfies the con- 
dition 77„j ^ 770 ^ 1 and spht the integral in H()8|) in two 
ones: E = Ei + E2 where 



2V2 



E2 — So 



2V2 



As it follows from Ht)7|) . 



dP 



1 (1-^^)' 



QY2 2P if 

and thus in the limit V ^ 



(We recall that the subindex "0" on the l.h.s. stands for 
V — and X — 0.) Finally, the definition of the effective 
mass H61(l . the explicit expression for the momentum Pq 
(|12|l . and the last formula H74|l yields the general expres- 
sion for the mass of the dark soliton: 



£0 



(75) 



In order to relate the mass of the soliton to the atomic 
mass m, we recall that —dE/dfi — N where N is the 
negative "total number of particles" associated with the 
soliton. Thus, tt^* = nis/N can be considered as the 
effective mass of a "solitonic" particle. In the limit V ^ 
the above quantities can be easily calculated to give 



dEi_^2V2 



1 

/ 



1 



a+1 



Vo 

d 1 



, ^ . V2 /- (l-r?^) dii 



£0 



V2 



a J drj \ Po[ri) / 77 0:770 
^ ^ 



(to obtain the last line we substituted the lower limit by 
zero, due to fast convergence of the integral, and inte- 
grated by parts). 

To calculate the derivative of E2 we take into account 
that ry is small over the whole range of integration. Thus 



no 



E2 



£0 



77 



2V2 



vo + l J J if - 



zdr] 



2V2 

£q , , 110 - £a 



2 V'n 



y/a+T a+1 rjo 

and in the limit V [due to ((HJ] , 
dE2 



dV^ 



-£o 



V2(a + 1) 



am 



(73) 



Sum of H72|) and (|73|) gives us the derivative we are look- 
ing for: 



dEp 
dV2 



— ^0 F^a ; ^a. 



1 d (1-77^) 

aj rjdri Pn{vi) 




-<i77. (74) 



N 



£0 



2a I (1 





drj 



and 



aFa 
771* = m . 



(76) 



(77) 



In Table HI we present examples of the effective mass 
for two relevant cases. From the provided values one can 
see that the effective mass of the soliton particle is bigger 
than the mass of a free particle: 771* > m. 



3. Frequency of oscillations of a dark soliton. 

To conclude this subsection we compute the frequen- 
cies of oscillations of dark solitons in a trap, what can 
be done using the relation Ht)2|) . To this end we notice 
that including the trap potential into the scheme devel- 
oped in the preceding subsection, can be done by simple 
change the chemical potential /j, by /j, — U{X). Thus for 
the parabolic trap we have 



dEo 



Uj"^ dEn 



2 ^^JL 

This leads us to the formula 

No. 



aFa 



(78) 



(79) 



B. Analysis based on the perturbation theory. 

In Ref. J^] it has been argued that the phenomenolog- 
ical approach formulated above can be justified from the 
viewpoint of the original GP equation with help of the 



perturbation theory for dark solitons |9|,|20j when motion 
occurs in a constant parabolic trap. The proof was based 
on a possibility of effective factorization of the solution 
on the constant background and the dark soliton solution 
moving against it. In the case of time dependent trap, 
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a 


So 


|7V| 






2 
1 


2^/31n^+/^£:o 
— ^° 


4^/31nl+/^ ^ 
V2 meg 


V3 + 2 In 

— m 

21ni±^ 

2m 


0.6572 a; 

LU 

71 



TABLE I: Characteristics (the zero-velocity energy Eq, the number of particles N , the effective mass m«, and the frequency of 
oscillations in the parabolic trap uja) of dark solitons with small velocities for different powers of the nonlinearity a. 



the background cannot be considered as a constant, and 
the theory requires revision. The goal of the present sub- 
section is to develop the modification of the perturbation 
theory and to obtain from it the exponent 7 which de- 
scribes change of the amplitude of oscillations of the GP 
dark soliton in a parabolic trap. 

To this end we start with the dimcnsionlcss form of 
the one-dimensional GP equation [the variable arc intro- 
duced in jnSl), see also Eq. 



takes the form 



roe. 



(85) 



We emphasize that the last equation is exact with no 
approximation made, so far. 

Before the analysis of (|85|l . let us consider Eqs. H83|l 
and (|84|l in more details. They can be reduced to a single 
equation for L: 



1 



(80) 



where v = ly (t) — h/ (2^/^CQm) lj (t) . We assume that 
CO (t) is a slow function of time, what is expressed by the 
adiabaticity condition ^ |^| 1- Accordingly, v (t) is 
also a slow function. Noliice that i/ (0) = i^o ^ 1 as a 
condition for the local density approximation. We look 
for a solution of l|80|) in a form of the ansatz (analogous 
of the well known lens transformation) 



1 



= '/'(C,T(r)), 



(81) 



where ^ is a funct ion o n time and on spatial coordinate 
given by ^ = C,/ ^/L(t)^ while f is a new temporal vari- 
able related to the old one by the equation ?> ~ \/L. 
The functions L{t) and /(r) are to be determined be- 
low. Substitution of (|HU into ((HOI yields 



i4>T + i'ii 



-^(L,-4/L)0-^(L,-8/L)e0c =0. (82) 

Let us now require the trap frequency of the new equation 
(i.e. term proportional to ^^0) to be constant, say 1/2i'q, 
and dissipative terms, i.e. the linear with respect to (j), 
to vanish. This gives us two equations: 



0, 



(86) 



Due to adiabaticity the first term in H8t)|l is small in com- 
parison with the other ones. Neglecting that term, we 
find in the leading order 



and f = -. 



(87) 



( 1 1^ 



Then, simple estimates give L^t 
and L vq/v ^ 1^ what justifies the approximation 
made. 

Next we introduce the notation R for the r.h.s. of 
(|85|l : R = —2ifL£^(f>^. Since / is small (because of the 
adiabaticity of the change of v) this term gives us an 
perturbation, which is a complementary to the perturba- 
tion introduced by a constant parabolic trap, provided 
1^0 is small, the case considered in detail in [9|. Due to 
their smallness, the effect of different perturbations on 
the dynamics of the soliton center is additive, allowing 
one to compute only the contribution of R to the dy- 
namical equation of the soliton center and add it to the 
equation describing soliton in a stationary potential ob- 
tained in |§| [see Eq. (32) there]. We skip description of 
tedious but straightforward calculations '2g\ and present 
only the final result: the equation for the soliton coor- 
dinate, in terms of the rescaled by L variables, is given 

by 



and 



2 , l.,2\^2 



fr + 4/- + -i. 



1 



(83) 



(84) 



The obtained equations will be supplied by the natural 
initial conditions /(O) = and i(0) = 1. Then Eq. ^ 



dr 2 



\vljx{T')dr 



VX. 



Next we differentiate the last equation with respect to 
T and eliminate the "dissipative" term by means of the 
substitution 



X = r(f)e''(^); 



(89) 
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where 



Then direct computation gives 



ln(^ 



(90) 



Having done this and restoring the original variables we 
arrive at the final formula 



1 , 



(91) 



Comparing this result with (|16|l and p8|l one can see 
a remarkable agreement. The perturbation theory, valid 
for relatively low densities of the condensate, and thus to 
the Gaussian background, corrects the law H18|l based on 
the phenomenological approach, by means of small shift 
(recall that V <^ 1 and thus 5 <C 1) toward larger expo- 
nent which in the TF limit is given by H16|l . Moreover, the 
perturbation theory introduces an explicit dependence of 
the exponent on the velocity (at this point it is relevant 
to recall also that the frequency itself does not depend 
on the soliton velocity) . Finally we mention that the ob- 
tained result corroborates with the numerical results on 
the dependence of the exponent 7 on the soliton velocity 
and on the density of the background depicted in Fig. ^ 
(d) and (h). 



J-xo c{x} 



(93) 



Thus for the small amplitude GP dark soliton in a 
parabolic trap we obtain to / v^, what coincides with the 
results known for relatively large velocity for the soliton. 
Eq. (|93|) alsogives cj^ = lj for a = 2, the result recently 
reported in |2a|. 

We emphasize however that presently there are no 
available results confirming validity of the law H93() for 
small amplitude NLS solitons. The main physical reason 
for this, mentioned in 0, is that in the vicinity of the 
turning points the density becomes small enough making 
the problem to be linear and thus not allowing solitonic 
propagation due to dominating dispersion. Mathemati- 
cally, the problem occurs due to divergence (see e.g. the 
second of equations (11) in Ref. [2^) of the small am- 
plitude expansion near the points where the condensate 
density, and thus the speed of sound, in the TF approxi- 
mation becomes zero (see also discussion of the failure of 
the small amplitude limit in p^l. 



VI. CONCLUSION 



C. A comment on small-amplitude solitons. 

Small amplitude dark solitons of the NLS equation 
with a polynomial nonlinearity of any power of the non- 
linearity are described by the KdV equation and they 
move with the sound velocity (or more precisely with 
a velocity slightly deviating from the sound velocity). 
While self consistent reduction of the 3D GP equation to 
one-dimensional KdV equation seems to be not possible 
for realistic condensates (as it is explained in Ref. [3|), 
the KdV being rather academic than practical allows one 
to predict some features of the underline GP equation. 

Moreover, one can easily argue, that the small- 
amplitude limit of a dark soltion in a parabolic trap is 
not available. Indeed, existence of a soliton in a trap 
implies smallness of the soliton width i compared to the 
trap width ^Jhj (muj). Using the expression for the width 
of a dark soliton, which is given by ^ , one immediately 
obtains the limitation huj/m < — u^. Thus the exis- 
tence of a trap does not allow the truth small amplitude 
limit, which would correspond to u — > c. 

Let now formally compute the half-period of oscilla- 
tions of a small-amplitude soliton in a parabolic trap. 
Under the half-period we understand the time necessary 
for a soliton to pass the distance between two turning 
points. To this end associate the velocity 



cix)=uj^-{xl~x') (92) 
where xq — 2/i/(ma;^), with the velocity of the soliton. 



In the paper we presented development of the theory 
suggested in the earlier publication [lOj, providing de- 
tailed description of the one-dimensional dynamics of a 
dark soliton in a Bose-Einstein condensate confined by an 
external potential. The theory is based on the local den- 
sity approximation and allows one to interpret the dark 
soliton as a hamiltonian particle. We addressed various 
generalizations of the theory including the nonlinearity 
of a general polynomial type as well as non-parabolic po- 
tential. We have obtained that the dependence of the 
amplitude of the soliton oscillations in a external trap 
depends on the adiabatically changing frequency through 
the scaling law Xq oc lu^'^ where the exponent 7 depends 
on the type of the nonlinearity and on the type of the 
confining potential. It turns out also that the frequency 
dependence of the amplitude of the oscillations depends 
also on the density of the condensate and on the initial 
velocity, even in the cases when the frequency itself is 
independent on the above quantities as in the case of 
the standard nonlinear Schrodinger dark solitons. Also 
the obtained scaling law in a general case appears to be 
very different from the corresponding law for the linear 
oscillator. 

We dedicated special attention to the cases of dark 
solitons within the framework of the Gross-Pitaevskii 
and quintic nonlinear Schrodinger models. We also have 
shown that in the limiting case of slow, and thus large- 
amplitude, solitons one can obtain the general explicit 
expressions for the effective mass of the dark soliton, con- 
sidered as a quasi-particle, and for the frequency of its 
oscillations in the external confining trap. 
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The results have been verified numerically, showing 
good agreement with theory, and were shown to be in 
agreement with outcomes of the direct perturbation the- 
ory for solitons. 
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In order to satisfy the constrain i^fe = at w = 0, we 
require Fr+i ^ Fr- From HA3|I and 1A4|I we obtain the 
following asymptotic relations 

Fk = O {Fr) = O {oj^'') , k<r 

which in their turn guarantee the smallness of the inte- 
grals 



dx oc uj'^'^ and 



d'^noix) 2r 

dx CX LU 



dx^ 



\x-X\<S 



\x-X\<S 



APPENDIX A: ESTIMATES FOR THE 
BACKGROUND 



when a; ^ 0. 



In the present appendix we provide estimates for the 
last two terms in (|37|l . For the sake of simplicity the con- 
sideration will be restricted to the case of a polynomial 
parabohc trap (p?^ . 

Let us consider the behavior of the function F{x) in 
the vicinity of the point x <^ Lq (we recall that Lq is an 
effective trap length). The background is obviously an 
even function of the trap what allows us to look for its 
solution in a form of the expansion 



APPENDIX B: ON THE LINK AMONG 
FORMULAS (inUl, (innil AND (IHHt 



In terms of the amplitude rj and the phase 9, both 
depending on ^ — 2Vt, Eq. 164|) can be rewritten in the 
form of a system 



F = l + J2FkC^ C^x^ 



(Al) 



fc=i 



More specifically we are looking for the coefficients Fk, 
all of which become zero in the homogeneous case when 
uj = 0. It follows directly from that 



= gn-Q Fi 



(A2) 



In the homogeneous condensate the chemical potential is 
given by /zq = griQ and thus there should be verified that 
Fi ^ mgriQ/h^ for u small enough. 

Next from one can obtain the recurrent formulas 

-nFk for fc < r (A3) 



-(2. + 2)(2. + l)^.,, = ^(-^ 



.^F, + _^2. (A4) 



2V'nc = "^vc^c + vScc 



2a 



(Bl) 
(B2) 



Multiplying (|B2p by rj, integrating with respect to 9 and 
using the boubdary conditions rj 1 and 9 — ^const as 
C — !■ ±oo, one obtains the link (|65|l . 

In order to obtain H66|) it is enough to substitute 0^ 
expressed in terms of rj through the relation (|65|1 in ljBl|) 
and multiply the result by rj. 



APPENDIX C: ADIABATIC INTEGRAL FOR 
THE GP SOLITON IN A POLYNOMIAL TRAP 



The adiabatic integral for the GP dark soliton is com- 
puted, using ((Tn|l and links and for a = 1, as 
follows 
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rgui 

^^l-l/(2r)^2+lA 

5^ 



arcsm 



(y2_^2)(2r-l)/(2r) y„2+^2 V 



where the constant Gr is given by 

22+l/(2r) .1 y(l+2/2) 



(l-y2) 



2\l-l/(2'-) 



y 



l + y2 
I 



arcsm 



dy. 



(CI) 



(C2) 



APPENDIX D: CALCULATION OF THE 
ENERGY (IH5t 

Starting with the definition H26|l . (|27|l . written as 
1 



and excluding with the help of (|12l) one obtains 
1 



E = Eq 



£; = 2En 



4Sn 



-(ry2"+2-l) + l-ry^ 



1 



a + 1 



„2q+2 



1 + 1 - 



drj 



iwr^2 



(Dl) Formula (|68() follows from the last equality. 
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